Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/cpp_libs/rwvector.lha/RWVector2.1/README

Text File | 1989-08-18 | 14KB | 339 lines

    The ROGUE WAVE Vector and Matrix Classes.
                Version 2

      Copyright (C) 1988, 1989.

      Dr. Thomas Keffer
      Rogue Wave Associates
      P.O. Box 85341
      Seattle WA 98145-1341

      Permission to use, copy, modify, and distribute this
      software and its documentation for any purpose and
      without fee is hereby granted, provided that the
      above copyright notice appear in all copies and that
      both that copyright notice and this permission notice
      appear in supporting documentation.

      This software is provided "as is" without any
      expressed or implied warranty.

*********************************************************************

    ******  A C K N O W L E D G E M E N T S  ******

I am indebted to Al Vermeulen (Univ. of Waterloo) for his insightful
feedback and debugging.  The critical idea of implementing all vectors
as slices was also his.

*********************************************************************

    ******  O V E R V I E W  ******

This is a set of vector and matrix classes.  Right now it includes all
standard arithmetic operations (+, *, /, -, --, etc.) for types
double, float, int, and complex, and most math functions (sin, cos,
sum, cumsum, etc.).  The Vector classes include three FFT "server"
classes:

1) DComplexFFTServer which can transform a complex vector
   into a complex vector (and back).

2) DoubleFFTServer which transforms a real vector into a complex
   conjugate even vector (and back).

3) DCosineFFTServer which performs sine and cosine transforms.  That
   is, it can transform either a real even or a real odd vector into a
   real even or odd vector (respectively).

The Matrix classes include an LU decomposition class (both double and
float) which can perform matrix inversions, solve systems of linear
equations, and calculate determinants.

The overall goal of the classes is to be as simple as possible while
not sacrificing (too much) in the way of speed.  Because it is easy to
generate absolutely huge header files and very long compilation times
when working with C++, simplicity is important.  Extendability is also
important.  Others are working on more exotic vector classes that
avoid the use of temporaries by building expression trees at runtime
and then optimizing them out.  No doubt these will be fast, but they
will also be very complicated, big, and hard to extend.

The classes here do not avoid using temporaries, but use reference
counting to minimize the time required for their construction.  They
should give you 90% of the speed of a more optimal package
(particulary when doing floating point intensive operations like
matrix inversion) with about a quarter of the size and compilation
times.  Example: the test program for LU decomposition (dludtest),
using g++, is only 90k bytes long (stripped).

The classes have been designed to be compatible with the NIH Vector
classes.  Indeed, they have the same names and many of the same
functions and syntax.  They implement "slices", but not "picks" and
"selections" (I have not found the last two very useful).  One
significant difference: a slice is taken with a slice() function,
rather than an overloaded operator().  I have found the latter
confusing.  Hence, use:
    a.slice(0,3,1)
instead of
    a(0,3,1)

The resulting syntax is very powerful.  For example:

    DGEMatrix a(4,4,0);           // A 4x4 matrix initialized to 0
    a.diagonal() = 1;             // Set the diagonal to 1
    a.diagonal(-1) = 2;           // Set the lower diagonal to 2
    a.row(2) = -1;                // Set row 2 to -1
    DGEMatrix ainv = inverse(a);  // Invert a

Note that many functions can be used as an lvalue.  This is also true
of vector slices.  Here's an example.  Use the Complex FFT Server to
verify Parseval's theorem at the Nyquist frequency.  It makes a good
example because it both uses slices and the FFT Server:

    DComplexFFTServer server;    // Construct the server

    // Make a test series, npts long, initialized to complex(1,1)
    DComplexVec Nyquist(npts, DComplex(1.0, 1.0));

    // Set every other element to complex(-1, -1):
    Nyquist.slice(1,npts/2,2) = DComplex(-1.0, -1.0));

    // Take its fourier transform:
    DComplexVec transform = server.fourier(Nyquist)/npts;

    // Calculate its original variance
    cout << variance(Nyquist) <<"\n";

    // Compare to the final, spectral variance:
    cout << spectralVariance(transform) <<"\n";

The use of slices can virtually eliminate the "for" loop over vector
and matrix elements.  This makes it easier to port the resulting code
to an array processory or other specialized hardware.  Or, a Basic
Linear Algebra (BLA) package can be used.

The test suite has many examples for using the classes.

*********************************************************************

    *******  Changes from Version 1  *******
        (not necessary to read right now)

The biggest changes from Version 1 is the addition of matrix classes
and the elimination of the Slice temporaries (e.g., "DoubleSlice").
Many small improvements were also made.

All vectors are now treated as slices.  This should be mostly
transparent to you.  This approach drastically reduces the number of
routines that must be written, the size of the code, and the number of
type conversions that have to be done at runtime.  It also makes type
conversions more reliable and obvious.  The one curious side effect
(as it turns out) is that you must pay attention to the very real
difference (in C++) between initialization and assignment:

    DoubleVec a, b, c;
    DoubleVec d = a;   // Note: d will be a REFERENCE to a
    d = b;             // Note: d will contain a COPY of b
    d.reference(c);    // Note: d will be a reference to d

Note that copy initialization is done by referencing the new variable
to the contents of the old variable, but assignment is done by COPYING
the contents over.  The function reference() is provided because a
variable can be initialized only once, subsequent assignments will
provide only copies.  It will simulate the effects of initialization.
The overall effect is that boths sides of an assignment MUST have the
same size.  Both variables need not have the same size when using
reference() (the object taking the reference will be resized to the
argument).

Why use a reference?  Speed.  They are used by the temporaries to
avoid unnecessary copying of the bulk of the vector.

*********************************************************************

    ******  I N S T A L L A T I O N  ******

The Classes can be compiled with either the GNU "g++" compiler
(version 1.35 or later) or the AT&T "cfront" compiler (Version 1.2, or
Version 2.0).  If the GNU compiler is used, then the class Complex, as
defined in the header file Complex.h distributed with the GNU g++
library is used for complex arithmetic.  If one of the AT&T compilers
is used, then class complex, as defined in complex.h is used.
However, the complex.h distributed with the V1.2 compiler cannot be
used to make vectors of complex: it has no unadorned constructor.  I
have included a revised version of complex.h in this distribution
which will automatically be used if you are using the 1.2 compiler.
The 2.0 version of complex has the necessary constructor.

0) Glance through the "Pitfalls" section below, just so you know
where the holes are.  Here's a language that's not even 5 years old,
and already there's no way to write a standard package!

1) Makefiles are provided for the GNU compiler (the default) and the
AT&T "cfront" compiler (Makefile.ATT).  Select the proper makefile.

2) AT&T compiler only:

   You might want to remove the #pragma statements in the files rw/*.h,
   if you don't want to look at warning errors.  There are two shell
   scripts for doing this: "depragmatize" to remove the #pragma statements,
   "pragmatize" to put them back in.

   Set the define near the top of Makefile.ATT to select either __ATT1__
   if you are using the 1.2 compiler, or __ATT2__ if you are using the
   version 2.0 compiler.

   Check the tail end of Makefile.ATT ("Conversions") to make sure
   that it is expecting the proper suffixes for your CC driver script.
   Mine accepts .C, puts out .C.o.  Yours might differ.

   You will have to set your environment variables as usual.

2) Make sure that the proper floating point options are selected for
your hardware in the appropriate Makefile.

3) As distributed, the GNU version 1.35.0 include file math.h does not
overload atan().  Version 1.35.1 does.  You might want to add this
entry.

4) The FFT functions depend on a fortran package from the National
Center for Atmospheric Research to do the actual FFT.  They are a lot
more complicated than some other routines, but they offer two
advantages: They let you do FFTs on vectors of arbitrary length (not
just powers of two), and they vectorize very nicely.  If you are
unable to compile the fortran routines (or don't want to!) then you
will be unable to run the full test suite.  Programs testdfft,
testcfft, and testcos require mathpack.  Or, you can supply your own
FFT routines.  A routine that transforms a complex vector into a
complex vector and back is all that is required.  All other transforms
(real to conjugate even, cosine transforms, etc.) sit on top of it.
These routines are put in a library called libmathpack.a.  If someone
has got a C++ (or C) function to do arbitrary length FFTs, I'd be
delighted to put it in.

5) The linear algebra stuff (LU decomposition, matrix inversion), etc.
are based on LINPACK, and also put in the library libmathpack.a (see
point 4 above).  Again, if you don't want to, or can't compile them,
then you will be unable to use the linear algebra routines.  Programs
dludtest and fludtest require these routines.  I decided to use
LINPACK because they are well know, robust, trusted, and have
assembler versions of their Basic Linear Algebra package.

6) cd to the src directory, then type make:

    cd src
    make

With the GNU compiler, you will get a slew of warning messages, due to
the inline nature of the GNU math library.  Here's a sample:

In function struct DoubleVec asin (const struct DoubleSlice &):
./../rw/DoubleVec.h:228: warning: argument passing of non-const * pointer from const *

7) To make the mathpack library (libmathpack.a), in the same src
directory, type

    make mathpack

8) To make the test suite, in the src directory, type

    make test

9) To verify it:

    make verify

I have tried to design the test suite so that variations at the
machine precision will not be apparent, but I may not have been
entirely successful.  In particular, the g++ and cfront formatting
differ slightly.

9) Finally, to install the package, check the Makefile macros LIBDIR,
and INCLDIR to see if the default installation locations suit you, then

    make install


**********************************************************************

    ******  P I T F A L L S  ******

* Most of the FFT routines *require* that
    sizeof(Complex) == 2 * sizeof(double)
in order to pass arguments correctly to the Fortran FFT routines.
This will be true provided:

    a) Structures can align on a double boundary (hard to imagine
    a machine where this isn't true);
    b) There are no other member objects in Complex besides the
    two doubles (i.e., real and imaginary);
    c) No base class in Complex;
    d) No virtual functions in Complex.
    e) Your compiler does not add any hidden data to the Complex
    class.  Neither the GNU nor the AT&T compilers do this.  Note,
    however, that this is NOT part of the language definition!

All of these conditions apply for both the GNU g++ Complex, and
Stroustrup's complex class and you shouldn't have any trouble.
Undefine the macro "COMPLEX_PACKS" in RComplex.h if one of these
conditions does not hold: Some of the routines have workarounds.
Others curl up and die.

* My version of the GNU loader (gcc-ld++) is unreliable in finding all
referenced externals.  For example, for a while it had trouble finding
"hypot.o" (which is in /usr/lib/libm.a).  If this happens to you, you
will have to extract the lost module and add it by hand to the TARGET
link list:

    ar x /usr/lib/libm.a hypot.o

(then change the TARGET entry in the Makefile to:)

${TARGET}:   ${TARGET}.o ${LIBFILE} ${MATHLIB}
    ${CCXX} ${CXXFLAGS} -o ${TARGET} ${TARGET}.o hypot.o ${LIBFILE} ${MATHLIB} -lm
                                                 ~~~~~~~

This bug has been reported.

* The file complex.h as distributed with the AT&T Version 1.2 compiler
cannot be used to make vectors of type complex.  This is because it
does not include an unadorned constructor complex::complex().  Hence,
I have modified the file and included it.  It is used only if the
Version 1 cfront compiler is used.  If the GNU compiler is used, then
the file Complex.h, as distributed with the GNU library, is used.  If
Version 2 of cfront is used then the version distributed with it is
used.

* Both compilers are touchy about inlines and type conversions.
Workarounds have been included.  In particular, the Version 1.2 cfront
compiler seems to have trouble with inlines nestled more than one
deep.

* I have not tried optimization flags with the routines.  I have no
idea if this would introduce compiler bugs.

* The AT&T version 1.2 compiler has trouble assigning a real to an
element of a complex vector.  Example:

    DComplexVec a(10);
    a(3) = 2.0;

will give a "bus error" (whatever that is).  The cure is to switch to
the GNU compiler.  No. Seriously.  Use a cast.  The last line should
be written as:

    a(3) = DComplex(2.0);


    ******  B U G   R E P O R T S  ******

Send 'em to:

 Dr. Thomas Keffer          | Internet: [e-mail]
 School of Oceanography     | BITNET:   keffer%sperm.ocean.washington.edu@UWAVM
 Univ. of Washington, WB-10 | uucp:     uw-beaver!sperm.ocean.washington.edu!keffer
 Seattle, WA 98195          | Telemail: T.KEFFER/OMNET
    [phone]